h4a_lm = lm_robust(attitudinal_outcome ~ h1_treat*islam, data = survey_df) %>% tidy()
h4b_lm = lm_robust(behavioral_outcome ~ h1_treat*islam, data = survey_df) %>% tidy()  

h6a_lm = lm_robust(attitudinal_outcome ~ h1_treat*past_attendance, data = survey_df) %>% tidy()
h6b_lm = lm_robust(behavioral_outcome ~ h1_treat*past_attendance, data = survey_df) %>% tidy()

h7a_lm = lm_robust(attitudinal_outcome ~ h1_treat*joko_supporter, data = survey_df) %>% tidy()
h7b_lm = lm_robust(behavioral_outcome ~ h1_treat*joko_supporter, data = survey_df) %>% tidy()

hte_mods = list(h4a_lm, h4b_lm, h6a_lm, h6b_lm, h7a_lm, h7b_lm)

hte_hypotheses = 
  bind_rows(hte_mods) %>%
  mutate(val = str_detect(term, "treat")) %>%
  filter(val == T) %>%
  mutate(val2 = str_detect(term, "\\:")) %>%
  filter(val2 == T) %>%
  mutate(outcome_lab = case_when(str_detect(outcome, "att") ~ "attitudinal compliance",
                                 str_detect(outcome, "beh") ~ "behavioral compliance")) %>%
  mutate(treatment_lab = case_when(str_detect(term, "islam") ~ "H4: Muslim \n x Any treatment",
                                   str_detect(term, "piety") ~ "H5: Pious \n x Any treatment",
                                   str_detect(term, "past") ~ "H6: Noncomplier \n x Any treatment",
                                   str_detect(term, "joko") ~ "H7: Jokowi Voter \n x Any treatment")) %>%
  mutate(treatment_lab = factor(treatment_lab, levels = c("H7: Jokowi Voter \n x Any treatment", "H6: Noncomplier \n x Any treatment",
                                                          "H5: Pious \n x Any treatment", "H4: Muslim \n x Any treatment")))
         
         


ggplot(data = hte_hypotheses) + 
  geom_hline(yintercept=0, lty=2) + 
  geom_pointrange(aes(x=treatment_lab, y=estimate, ymin = (estimate - 1.96*std.error), ymax=(estimate + 1.96*std.error))) + 
  geom_pointrange(aes(x=treatment_lab, y=estimate, ymin = (estimate - 1.65*std.error), ymax=(estimate + 1.65*std.error)), fatten=2, size=1.1) + 
  labs(x='', y = '') + 
  ylab(expression("Conditional Average Treatment Effect (" ~ beta[3] ~ ")")) +
  facet_wrap(outcome_lab~.) + 
  coord_flip() + 
  theme_bw() +
  theme(legend.title=element_blank(), text=element_text(size=20), 
        axis.title.y=element_blank()) +
  theme(text=element_text(size=12, family="Times")) + 
  theme(panel.grid.minor = element_blank(), panel.grid.major.x = element_blank(), panel.grid.major.y = element_line(size = .2)) + 
  theme(strip.background =element_rect(fill="white"))


ggsave("./outputs/main_figure_2.pdf", width = 6, height = 3)

rm(list=setdiff(ls(), "survey_df"))
